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Abstract 

A lubrication model describes the dynamics of a thin layer of fluid 
spreading over a solid substrate. But to make forecasts we need to 
supply correct initial conditions to the model. Remarkably, the initial 
fluid thickness is not the correct initial thickness for the lubrication 



model. Theory recently developed in [12, 14] provides the correct 
projection of initial conditions onto a model of a dynamical system. 
The correct projection is determined by requiring that the model's 
solution exponentially quickly approaches that of the actual fluid dy- 
namics. For lubrication we show that although the initial free surface 
shape contributes the most to the model's initial conditions, the ini- 
tial velocity field is also an influence. The projection also gives a 
rationale for incorporating miscellaneous small forcing effects into the 
lubrication model; gravitational forcing is given as one example. 
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1 Introduction 

The flows of thin films of fluids are encountered in many engineering and 
biological applications. They include: the flow of rainwater on a road or 



windscreen or other draining problems 0; paint and coating flows |T7|, [1 
the flow of many protective biological fluids and other coating, painting 
and drying processes [|K], |], ^(J, e.g.]. The fluid film thickness and the average 



fluid flux are the main characteristics of interest in these applications. The 
fine details of the actual local velocity and pressure fields typically are of less 
practical importance. For this reason the various approximations have been 
constructed over the past decades 0, [U], [18], [13|, [7|, e.g.] to model the evo- 



lution of the fluid flow in various geometries |16[ and parameter regimes ||. 
We consider herein the basic nondimensional lubrication model for surface 
tension dominated flows, 



Vt 



\dx (v 3 Vxxx) , (1) 



where r](x, t) is the thickness of the fluid film spreading over a solid substrate 
(at y = 0). Centre manifold theory [|TJ provides a generic and systematic pro- 
cedure for deriving such models for a wide variety of fluid flows. Recently, 
Roberts et al. [[15], [LJJ showed how this well established lubrication model of 
thin film flow is rigorously derived from the governing Navier- Stokes equa- 
tions (outlined in Section |2|) using a computer algebra implementation of 
centre manifold techniques. The model is derived under the assumption that 
longitudinal derivatives, d/dx, are small — the slowly varying assumption — 
as used extensively in creating models of shear dispersion, || [L9|, e.g.]. The 
centre manifold based algorithm provides a straightforward derivation of the 
model up to arbitrarily high order [13], but in this work we limit ourselves 



to consideration of the above leading order model. 

It is frequently believed that the initial conditions for the models of long- 
term evolution are not highly important because the asymptotic state does 
not depend on a character of transient processes in the system. It was shown 
in Hl2|, M that initial conditions can have a long-lasting influence on forecasts. 



This seems especially true for models of spatio-temporal dynamics such as 
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the lubrication model. Similarly, in the case of dispersion in a channel |L2| or 
pipe @ the long term location and spread of a pollutant does depend upon 
the details of the initial release of the pollutant. Thus the initial conditions 
for a model must be chosen carefully to ensure the long term fidelity between 
the model and the physical flow. It will be shown in Section ||] that the initial 
form of the free surface for the lubrication model should be different from the 
initial thickness of the physical fluid. In particular, if the fluid initially has 
thickness f]o(x) but zero velocity and pressure then the lubrication model ([!]) 
should be solved with initial condition that r](x, 0) = h (x) where 

h ^r]o + VoVoxx ■ (2) 

In general, the initial condition ho for the model is the non-trivial func- 
tion of initial velocities and pressure distributions given by fl3"I|)-(|3l)D . The 
argument for these initial conditions is based upon the dynamics near the 
low-dimensional centre manifold, that is, upon the physics of the approach 
to the lubrication model. The general arguments, developed in |12| and 



recently refined in [14], are based upon the geometric picture provided by 
centre manifold theory. The principle aim of this paper is to apply this gen- 
eral framework to the considerable complications of the infinite dimensional 
dynamics of thin film fluid flow and so derive (0) and its generalisations. 
This is the first time that correct initial conditions have been obtained for 
lubrication theory. 



An interesting aspect of the general analysis developed in |12j and [14 
is that the projection of initial conditions also gives a rationale for treating 
small forcing of the dynamics. This connection was more fully explored by 
Cox & Roberts who discussed the effects upon the centre manifold and 
the evolution thereon for time dependent forcing. In Section ^ we apply 
the projection to a gravitational forcing of the thin fluid layer to verify the 
veracity of the classic model 

Vt-~d x r] 3 (r] xxx + Bsm9-Bcos9r] x ) , (3) 

where B is the nondimensional magnitude of gravity and 9 is the downwards 
angle of the substrate. This is derived here as just one example of a very 
general result that applies to all small forcing effects upon the fluid flow. 



2 The lubrication model of fluid flow 



We consider a two-dimensional flow of thin film of Newtonian fluid along a 
flat horizontal substrate. The free surface is given by y = rj(x, t), where x and 
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y are horizontal and vertical coordinates respectively. The flow, with velocity 
q = (u, v) and pressure p, is governed by the incompressible Navier-Stokes 
equations 

q t + q ■ Vq = —Vp + z/V 2 q , (4) 
P 

supplemented by the continuity equation 

V-q = 0, (5) 
non-slip boundary conditions on the bottom 

q = on y = , (6) 
and tangential stress and normal stress conditions on the free surface 

2r]x(vy - u x ) + (l - r£*j (u y + v x ) = on y = r] , (7) 



(l + rj^j p = 2/i v y + rf x u x - r] x (u y + v x ) 



on y = r] , (8) 



1 + 



respectively, as discussed in detail by Roberts [|T^|. We close the problem 
with the kinematic condition relating the velocity of the fluid on the surface 
to the evolution of the free surface: 

f] t = v -urj x on y = 7] . (9) 

In the above equations p is the fluid density, u is the kinematic viscosity, and 
a is the coefficient of the surface tension. The fluid film is assumed to be so 
thin that the gravity force in the momentum equations can be neglected (see 
the discussion in Hl5| ) at least initially. 



We non-dimensionalise the governing equations using a typical film thick- 
ness if as a reference length, reference time [iHjo (where \x is the dynamic 
viscosity of the fluid), reference speed cx//i, and reference pressure a/H. On 
this small scale, fluid viscosity is strong and the fluid layer is of very large ex- 
tent laterally. The non-dimensional Navier-Stokes and continuity equations 
then become 



?e(q t + q-Vq) = -Vp + V 2 q, 
V-q = 0, 



(10) 
(11) 
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where 1Z = paH/fi 2 is a Reynolds number. They are complemented by the 
normal stress condition 

f 1 + t£) P = 2 K + vlu x - r] x (u y + v x )] Vxx on y = rj , (12) 

V ' L J ^ + V 2 X 

and equations @, (0) and (§) which remain symbolically unchanged under 
the nondimensionalisation. 

In such a thin layer of fluid, the infinite number of horizontal shear modes 
decay exponentially quickly through viscous dissipation acting across the 
thin film. Thus in the long term, the dynamics are driven by surface tension 
trying to flatten surface curvature. Centre manifold theory JTJ is used in 
such circumstances to systematically derive the low- dimensional model of 
the long term evolution, here the lubrication model (|I|) for the fluid layer's 
thickness 77, see fL5| , §3] for more introductory detail. The approximate form 
of the lubrication model for such a flow is obtained as a formal expansion 
in orders of the x-derivatives under the assumption that these derivatives 
are small. Although the model can be developed to an arbitrary order of 
spatial derivatives using the iterative computer algebra algorithm suggested 
in Jl5]| , the expressions for higher order approximations are very involved and 



thus we present here only the lowest order model. To errors of fifth-order 
in d x and parameterized by the free surface thickness 77, the centre manifold 
v = (u(r]),v(r]),p(r]),7]) is given by 

u « (yj) - ]-y 2 ) r] xxx , (13) 



2 

Vo 2 



1 y 2 VxVxxx + (lv 3 -lv 2 v)d x v, (14) 



p ~ -Vxx + ^vlvxx -(v + y) VxVxxx - {^rj 2 + yv- \y 2 ^j dtv , (is) 

where 77 evolves according to ([I]). Observe that up to this order the model 
does not depend on the Reynolds number — fluid inertia is negligible. The 
lubrication model (0) is the basic model for the dynamics of thin fluid films. 



3 Project the initial conditions 

In order to use the lubrication model ([j]) to make forecasts, it should be sup- 
plemented with initial conditions. Roberts |12| has shown that determining 
the correct initial conditions is a nontrivial problem. Remarkably, in general 
the initial value of 77 for model ([I]) differs from the initial fluid thickness for 
the physical problem (p|)-(|T2|). To distinguish between the two, we denote 
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the initial fluid thickness by r] and use h to denote the initial conditions 
for model (|l|) of the fluid's evolution. The main task of this paper is to show 
how to determine ho as a function of the initial fluid state. 

To define the proper initial conditions for model ([I]) we follow the pro- 



cedure outlined by Roberts |14| and examine the dynamics in the vicinity 
of the centre manifold. We start by linearizing the governing equations 
about the centre manifold (|i3l)-(|I"5|) by writing the fluid variables as the 
sum (u, v, p, rf) + (V, v ', p', rj') where primed quantities are the assumed small 
displacement from the centre manifold, and so their products are neglected. 
The resulting system is 

ft (q* + q' • Vq + q • Vq') + Vp' - V 2 q' = 0, (16) 

Vq = 0, (17) 

with the boundary conditions at y = rj 

r]' t -v' + u'r) x + ur]' x + u y r) x r)' - v y r\ = , (18) 

2 Vx (v y ~ Ux) + 2% (y' y - u x ) + (l - vTj {u'y + <) - ^VxVx + v x ) 

+2r] x (v yy - u xy )r]' + (l - r£) (u yy + v xy )r]' = , (19) 

2\ „/ i o„ „' „ i ft i „2\ „ „/ , Vxx VxxVxVx 



(l + r] 2 x ) p' + 2r] x T]' x p + (l + rf x ) p y rj + 
2 



1+7/2 (1 + ^)3/2 

y + Vlu' x + 2VxV x Ux - Vx (u'y + <) - Vx ( u y + v x)] (20) 



+2 



'^Vy Vx^Xy Vxi^Uyy + U x y / 



rf 



and the homogeneous boundary conditions for the velocity q' = at y — 
0. The above equations describe the dynamics of the fluid near the centre 
manifold ©-(ll). 

Typically, the initial conditions u = (u (x, y), v (x, y),p (x, y), 7/0(2)) 
for the original fluid layer equations (||)-(|1^) do not belong to the low- 
dimensional centre manifold v given by (fl3| ) -(|T5l ) . Thus they cannot be 
used directly as a starting point for model ([[]). As shown in |12| and |L4[] the 
proper model initial condition is the projection v = (u(h Q ),v(h ),p(h ), h ) 
from Uo to the centre manifold along the isochron — in the state space an 
isochron is a surface of all the initial states which have the same long-term 
dynamics on the centre manifold (up to an exponentially small error). Con- 
sequently, the model initial conditions are determined to satisfy 



(z, u - v ) = , 



(21) 



}3: Project the initial conditions 



where z = {u\v\p\ry) is a vector orthogonal to the direction of projection 
(the dagger is used to denote field quantities in the adjoint space). Here the 
inner product is defined for four component vector fields 

a= (a 1 (x, y, t), a 2 (x, y, t), a 3 (x, y, t), a 4 (x, t)) and 
b = (b 1 (x,y,t),b 2 (x,y,t),b 3 (x,y,t),b A (x,t)) 



as 



(a,b) 



OO ft] 



oo JO 



a L b L + a 2 b 2 + a 6 b A ) dydx + / a 4 6 4 dx . 



(22) 



According to the arguments developed in |T2| and refined in fl4|l , the defining 
vector of the projection, z, satisfies the dual equation 



Vz = (Vz, e)z , 
where e is the local tangent vector to the centre manifold 



+ 0(d: 



(23) 





u 




" " 


d 


V 







dr] 


V 




-d 2 




. V . 




1 



(24) 



and the dual operator V is obtained from equations adjoint to (p^)-(pO|) with 
respect to the inner product defined as 



(a, b» 



(a, b) dr 



(25) 



Higher order derivative terms in ( p4|) can be easily computed but, as will 
be shown later, the given second order truncation will suffice for finding the 
initial conditions to the first few orders. Note that the local tangent to the 
centre manifold is a vector operator rather than just a vector function as oc- 



curs in the finite dimensional cases discussed in |T^]. Being introduced into 
the inner product (0), it acts on the other vector involved before the integra- 
tion is performed. Using the above inner products, the adjoint expressions 
of (|I~6|)-(|2"0|) leading to the dual operator T> are: 



Vz 



1Z (u\ — v)u x — v^v x + uul + vv)^j + pi + u xx + ujy 
TZ (y\ - u^u y - vh y + uv\ + vvl) +pl + v xx + v\ 

ul + v\ 

Vt + ~ P ] + v^Vxx + ul + 2v ] x r] x -vl on y = r] 



yy 
yy 



(26) 



The adjoint velocities satisfy homogeneous boundary conditions q' = at 
y = 0. The adjoint boundary conditions for the velocities at y = r] are 
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represented by expressions too long to be given here. Instead in the Appendix 
we include computer algebra code written in reduce for obtaining them and 
the dual operator D. Periodic boundary conditions at x — ±00 are used in 
the derivation. 

Given the above dual T>, system ( ^3|) is solved asymptotically assuming 
that it is possible to neglect higher order derivatives with respect to x. The 
treatment of d x as small is equivalent to the assumption of slow variation in 
x. The iterative algorithm is quite similar to the one described in |T3[ and 
used to derive the centre manifold model (|13|)- (|15|) & (jl]). Thus here we just 
make a few notes on the specifics of the first few iterations. In essence the 
procedure is as follows: we start by solving the equations neglecting all x 
derivatives and then in further iterations compute the corrections associated 
with these derivatives of functions found at previous iterations. Owing to 
the special form of the vector e, the right-hand side of (|23| ) remains zero 
during the first few iterations required to obtain the leading order (in d x ) 
expressions for z. It is easier to first look for the solution z in the functional 









' 




«t 







z = 


pi 




7y 




77T 




T) 



s.t. rji 



+ 



0. 



+ 



vl (vv - \y 2 






Vxx 



Vxx 



(I?/ 3 - \y 2 v) - \vlvxy 2 



(yv - \y 2 ) + vlvxy + \d x (r^ 2 ) 



(27) 



(2f 



From the structure of the successive corrections to the solution for z we 
deduce that the initial conditions for the model are influenced the most by 
the initial form of the fluid surface. This is the expected result for such a 
surface tension dominated flow. The initial horizontal velocity field has a 
secondary effect on the flow (corresponding terms in ( p7[) appear only in the 
second iteration) primarily as a response to the horizontal pressure gradient 
induced by the surface curvature. The vertical motion is even less important 
since it is severely restricted by the small thickness of the fluid layer. 

An additional condition must be exploited to determine as an asymp- 
totic expansion in d x . It is done using the normalisation condition (z,e) = 1 
(see WW) which upon (El) leads to 



Pxx ' 



dy) dx + O (dl) 



(29) 
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At the leading order we then obtain 



/oo 
rf dx = 1 . 
-oo 



(30) 



This condition does not provide any unique solution for but rather a 
continuum of linearly independent localised functions. Without any loss of 
generality we choose the linearly independent solutions 



rf(x; a;*) = 8{x - x*) + O (d^j 



(31) 



where x* is an arbitrary point and 5 is the Dirac delta function. At the next 
iteration we require the second order contribution to to vanish. This 
results in 



rf(x; x*) = 5(x - x*) + 5"(x - x*)r](x, t) + o (d^j . 



(32) 



This expression for is used in (^) to determine z, the defining vector of 
the proper projection onto the centre manifold. 

Lastly, we use z in to project an initial condition. Requiring the inte- 



grand in equation fl2~T| ) to vanish and taking into account 
we obtain 



and (P|) 



T] - h + p + d x 



u y 



| -h 



hox (po(h + y) - \v y 2 



h n h 



O'l-Oxx 



(33) 



+dl 



hoVo -hl + po (^*- + h (y + 1)) - f (h y 2 - 







where the notation / = J m f dy is introduced. This equation determines ho 
and can be solved iteratively as well. The first three iterations produce 



ho 



hoo + hi + h, 



iQ2 



(34) 



where 
hoo 
hoi 



Vo+Po 



II 



h 



02 



-d x [h 01 Uoy] + h 00 h oxx - d x 
+dl po 



(35) 



''OOx 



Po(h t 



00 



y) - 2 v ov 2 



'h\o-y 2 



v 



h 00 y J - — I h 00 y 2 - — 



Note some specific cases of interest. 
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• Parallel shear flow v = p = 0, r] Q = 1 and u = u (y). Accord- 
ing to the linear stability analysis, due to viscous dissipation such 
a flow approaches the motionless state exponentially quickly Thus 
the centre manifold model, which disregards exponentially fast tran- 
sients, must give rise to a stationary solution. In this case the initial 
conditions ([34]) for the model and the model solution itself are just 
v = v(t) = (0, 0, 0, 1) and do indeed correspond to the motionless 
uniform fluid film. 

• Initially stationary fluid layer (w = ^0 = 0) with uniform pressure 
equal to the atmospheric one (p = 0) and curved free surface r/ ^ 
const. In this case ho « rj (l + rjo xx ) and the initial conditions for 
the model coincide with the initial film thickness only to leading or- 
der. Higher order terms tend to smooth out the initial distribution of 
the model film thickness flattening "hills" and "valleys" . This can be 
interpreted in the following way. The physical fluid which is initially 
motionless requires time for acceleration i.e. time to approach the centre 
manifold ([T3|) — (p~5|) in which velocities are generally non-zero. Dissipa- 
tion during this transient acceleration leads to a decrease in the energy 
of the system. Since the energy of the system up to leading order in 
d x is just the potential energy associated with the surface tension and 
is proportional to the surface curvature the initial condition models 
the energy loss by levelling out the free surface in comparison with the 
original distribution. 

• Nonzero initial average pressure. It leads to a change in the initial 
model fluid film thickness when compared with that of the original 
problem. In particular, locally positive initial pressure corresponds to 
locally thicker fluid film in the model. This is intuitively expected 
since the increased pressure inside the fluid layer (imagine an under- 
water explosion) acts against the local surface tension and leads to the 
appearance of the local "hill" on the film surface. 

Finally we note that the initial condition for the model is most sensitive 
to the fluid film thickness and the local pressure whereas the initial velocity 
field has just secondary effect on the long term film dynamics. This is not a 



surprise since, as noted in |15| , the considered flow is essentially the creeping 
one and inertia effects are less important than the influence of the surface 
tension or, equivalently, of the surface curvature. 
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4 Gravitational forcing as an afterthought 



In the previous sections the influence of gravity on the thin film flow is 
neglected. Here we demonstrate how such a forcing may be added into the 
model using the projection derived for initial conditions as argued in general 
in [14]. The general technique may be used to include physical processes into 
the lubrication model after developing the model. The result given here for 
gravity is just one specific example. 

The correction to the model accounting for the gravity can be obtained 
by iterative solution of the Navier-Stokes equations ( |10D where terms f = 
B(gi,g2) responsible for gravity are introduced in the right-hand side [ IE , 
e.g.]. Here B = p\g\H 2 /a (assumed finite but small) is the Bond number 
flT5| and gi and #2 are components of the non-dimensional gravity vector in x- 
and y-directions, respectively. Alternatively, considering gravity as a specific 



example of a forcing which by arguments in |14j can be directly projected 
onto the model ([I]). Geometrically, the centre manifold obtained for the 
dynamics without forcing is deformed slightly when forcing is applied such 
that each point of the unforced centre manifold is shifted along the isochron 
passing through the original location of this point as discussed by Cox & 
Roberts ||. According to the dynamics of the free surface subject to 
forcing is described by the modified model (P 



Vt 



1 ( 3 
~7^9x \ V Vxxx 



where q is the projection of the forcing f of the fluid, namely 

g=(z,f). 



(36) 



(37) 



Typically, gravity is uniform in thin film flow applications and then (gi, #2) = 
(sin 8, — cos#), where 9 is the downwards angle between a flat substrate and 
the horizontal. Then upon using (|27D and (|3l|) , (|37|) immediately leads to 



Q 



B 



1 

+ 2 92 



Ul <t 1 '/ -y'J 



y 



-y 

3 y 



S'xVx 



y 2 > dy dx 



(3? 



-B 



giv 2 Vx + -gidx 



which is identical to the correction obtained from directly modelling the 
forced equations (P~0| ) through assuming small Bond number \TE\. The theory 



of initial conditions recently refined in |TJ] also provides an elegant way of 
modifying the model to incorporate forcing. 
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5 Conclusions 

The proper initial conditions for the lubrication model of flow of thin film is 
derived using the projection of the initial conditions for the original problem 
onto the centre manifold representing the lubrication model. The obtained 
results are easily generalised to the case of isotropic three-dimensional thin 
film flow. Then the two-dimensional lubrication model 



should be solved with initial conditions given up to the first order by 



Here qo = (uq,wq) is the initial horizontal velocity field for the original 
problem and V is a two-dimensional operator in xz-pl&ne. 

Acknowledgement This work was supported by a grant from the Aus- 
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A Computer algebra derives the dual 

reduce^ code for determining the dual operator T> along with its associated 
boundary conditions: 

1 °/ Seek to find adjoint of thin fluid film equations. 

2 °/ Linearise about known centre manifold vu,vv,vp,vh. 

3 °/ a denotes adjoint quantities, b - physical ones. 

4 factor b; on div; off allfac; °/ improves the output text 

5 operator b,a; 

6 depend b,x,y,t$ depend a,x,y,t$ depend vh,x,t$ 

7 depend vu,x,y,t$ depend vv,x,y,t$ depend vp,x,y,t$ 

8 % look at linearisation about centre manifold dynamics 

9 let del~2=>0$ % gets rid of all nonlinear terms 

10 qu:=vu+del*b(u)$ qv:=vv+del*b(v)$ 

11 pp:=vp+del*b(p)$ hh:=vh+del*b(h)$ 

12 % some shorthands 

13 depend rsqrt,x,t$ % short for 1/sqrt (l+h_x~2) 

*At the time of writing, information about reduce was available from Anthony 
C. Hearn, RAND, Santa Monica, CA 90407-2138, USA. E-mail: reduce@rand.org 




(39) 




(40) 
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14 depend rmin,x,t$ % short for l/(l-h_x~2) 

15 let { df (vh,x)~2*rsqrt~2 => l-rsqrt~2 

16 , df (vh,x) ~2*rmin => rmin-1 

17 , (df (vh,x)~2+l) => l/rsqrt~2 

18 , df (rsqrt , ~zz) => -rsqrt~3*df (vh,x)*df (vh,x,zz) 

19 , df (rmin, ~zz) => 2*rmin~2*df (vh,x) *df (vh,x,zz) 

20 , rmin*rsqrt~2 => (rmin+rsqrt~2)/2 }$ 

21 % physical equations 

22 umom:=re*(df (qu,t)+qu*df (qu,x)+qv*df (qu,y) )+df (pp,x) 

23 -df (qu,x,x)-df (qu,y,y)$ 

24 vmom:=re*(df (qv,t)+qu*df (qv,x)+qv*df (qv,y) )+df (pp,y) 

25 -df (qv,x,x)-df (qv,y,y)$ 

26 cty:=df (qu,x)+df (qv,y)$ 

27 f kin : =df (hh , t ) -qv+qu*df (hh , x) $ 

28 f tan : =2*df (hh , x) * (df (qv , y) -df (qu , x) ) 

29 +(l-df (hh,x)~2)*(df (qu,y)+df (qv,x))$ 

30 fnor:=(l+df (hh,x)~2)*pp 

31 -2* (df (qv , y) +df (hh , x) ~2*df (qu , x) 

32 -df (hh , x) * (df (qu , y) +df (qv , x) ) ) 

33 +df (hh , x , x) /sqrt ( 1+df (hh , x) ~2) $ 



34 % linearization of boundary and kinematic conditions 

35 let{df (vv,y)=>-df (vu.x) ,df (b(h) ,y) =>0}$ 

36 f kin : =df (f kin+del*b (h) *df (f kin , y) , del) $ 

37 rll:={df (b(v) ,y)=>-df (b(u) ,x)}$ 

38 let rll$ 

39 ftan:=rmin*df (ftan+del*b(h)*df (ftan,y) ,del)$ 

40 fnor:=rsqrt"2*sub(del=0,df (f nor+del*b(h) *df (fnor.y) ,del)) ; 

41 clearrules rll$ 

42 let {df (vh,t)=>vv-vu*df (vh,x)}$ 

43 °/ linearized equations 

44 umom:=df (umom,del)$ vmom:=df (vmom,del)$ cty:=df (cty,del)$ 

45 °/ innerproduct form and integrate 

46 on list$ 

47 operator iii$ linear iii$ 

48 depend x,xyt$ depend y,xyt$ depend t,xyt$ 

49 let {iii(~aa*df (b(~bb) ,y) ,xyt) => 



50 -df (aa,y)*b(bb)+fs*aa*b(~bb) , 

51 iii(~aa*df (b(~bb) ,x) ,xyt) => 

52 -df (aa,x)*b(bb)-fs*df (vh,x)*aa*b(~bb) , 

53 iii(~aa*df (b(~bb) ,t) ,xyt) => 

54 -df (aa,t)*b(bb)-fs*df (vh,t)*aa*b(~bb) , 
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55 iii(~aa*df (b(~bb) ,y,2) ,xyt)=> 

56 -iii(df (aa,y)*df (b(bb) ,y) ,xyt)+f s*aa*df (b(~bb) ,y) , 

57 iii(~aa*df (b(~bb) ,x,2) ,xyt)=> 

58 -iii(df (aa,x)*df (b(bb) ,x) ,xyt) 

59 -f s*df (vh,x)*aa*df (b(~bb) ,x) , 

60 iii(a(~aa)*b(~bb)*~cc,xyt) => a(aa)*b(bb)*cc}$ 

61 iadj :=-(iii(a(u)*umom+a(v)*vmom+a(p)*cty,xyt) 

62 +f s*a(h)*fkin)$ 

63 °/ extract adjoint PDEs 

64 umom:=df (sub (fs=0, iadj) ,b(u))$ 

65 vmom: =df (sub(f s=0, iadj ) ,b(v) )$ 

66 cty:=df (sub(f s=0,iadj) ,b(p))$ 

67 °/ extract the adjoint FS boundary conditions 

68 let {df (b(v) ,y)=>-df (b(u) ,x)}$ 

69 °/ df(u,y)=buy on the surface 

70 depend fl,x,t$ depend f2,x,t$ depend f3,x,t$ 

71 depend f4,x,t$ depend f5,x,t$ depend f6,x,t$ 

72 depend f7,x,t$ depend f8,x,t$ depend f9,x,t$ 

73 buy:=-df (b(v) ,x)+f l*df (b(u) ,x)+f 2*df (b(h) ,x)+f3*b(h)$ 

74 bp : =f 4*df (b (u) , x) +f 5* (df (b (v) , x) +buy) +f 6*b (h) +f 7*df (b (h) , x) 

75 +f8*df (b(h) ,x,2)$ 

76 iadj :=sub(b(p)=bp,sub(df (b(u) ,y)=buy,df (iadj ,f s)))$ 

77 operator ii$ linear ii$ 

78 depend x,xt$ depend t,xt$ 

79 let {ii(~aa*df (b(~bb) ,x) ,xt) =>-df (aa,x) *b(bb) , 

80 ii(~aa*df (b(~bb) ,y) ,xt) => aa*df (b(bb) ,y) , 

81 ii(~aa*df (b(~bb) ,t) ,xt) =>-df (aa,t)*b(bb) , 

82 ii(~aa*df (b(~bb) ,x,2) ,xt) =>df (aa,x,2)*b(bb) , 

83 ii(~aa*b(~bb) ,xt) => aa*b(bb)}$ 

84 iadj :=ii (iadj ,xt)$ 

85 factor a$ 

86 bl:=df (iadj,b(h))$ b2 : =df (iadj ,b(u) ) $ b3:=df (iadj ,b(v))$ 

87 % define coefficients entering the definition of buy 

88 f l:=-coeffn(ftan,df (b(u) ,x) ,1)$ 

89 f2:=-coeffn(ftan,df (b(h) ,x) ,1)$ 

90 f3:=-coeffn(ftan,b(h) ,1)$ 

91 f4:=-coeffn(fnor,df (b(u) ,x) ,1)$ 

92 f5:=-coeffn(fnor,df (b(v) ,x) ,1)$ 

93 f6:=-coeffn(fnor,b(h) ,1)$ 

94 f7:=-coeffn(fnor,df (b(h) ,x) ,1)$ 

95 f8:=-coeffn(fnor,df (b(h) ,x,2) ,1)$ 
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5 % output for the resulting equations 
7 off nat; out "adj. red"; 
3 umom; vmom; cty; 
9 bl; b2; b3; 

3 write "end"; shut "adj. red"; on nat; 

1 end$ 

2 

3 

Below we comment on the listed REDUCE program: 

1. Preliminaries. 

• £ 10-12 describes states slightly displaced from the centre mani- 
fold. 

• £ 13-20 defines short hands and their properties for the expressions 
entering the free surface boundary conditions and appropriate al- 
gebraic and differential rules. 

• £ 22-33 expresses the physical fluid equations and their boundary 
conditions. 

2. The linearisation about the centre manifold. 

• £ 35 makes use of the continuity equation to get rid of v y and 
affirms that the free surface form rf does not depend on the vertical 
coordinate y. 

• £ 37, 38 state that v' = —u' x in order to simplify the boundary 
conditions. This definition has to be local to allow for the deriva- 
tion of the adjoint continuity equation later. Thus £ 41 clears this 
rule after it is used here. 

• £ 36, 39 and 40 extracts the linearized kinematic and boundary 
conditions taking into account the variation in y of the free surface 
itself. 

• £ 44 extracts the linearized momentum and continuity equations. 

3. Determination of the adjoint equations. 

• £ 47 introduces the operator iii which obtains the adjoints to the 
differential sub-operators entering linearized equations (p!6|)- ([20|) 
through the integration by parts rules listed in £ 49-60. Note 
that the volumetric integrals in (|22"D after integration by parts 
contribute to the surface integral because the adjoint functions 
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and their derivatives generally are not zero on the free surface 
which is a function of x and t. In addition, note that the rule in 
£ 53-54 uses the rule previously defined in £ 42. 

• £ 61-62 forms the inner product ( p2|) taken with a negative sign 
for further convenience. 

• Finally, £ 64-66 extracts the adjoint equations which result di- 
rectly in the expression for the dual to be solved to obtain the 
model initial condition generating functions. 

4. Determination of the adjoint kinematic and boundary conditions. 

This is done in three steps: 

• Firstly, the y derivatives of the unknown functions must be elim- 
inated from the expressions for the adjoint boundary conditions 
since they remain undefined under the surface integration. This 
is done by making use of the continuity equation to eliminate 
v' (£ 68-76) and the tangential stress boundary condition {£ 73 
with so far not defined coefficients fi) to eliminate u' y . Secondly, 
p'(x,r] + i]',t) is eliminated through the normal stress boundary 
condition {£ 74-75). 

• The surface operator ii is introduced in £ 77, which specifies the 
integration by parts rules {£ 79-83) along the free surface. The 
adjoint boundary conditions are obtained by acting with the op- 
erator ii on the redefined in £ 76 inner product iadj {£ 84-86). 

• Finally, the undetermined coefficients fi are determined in £ 88-95 
and the final output is written in the separate file {£ 97-100). 
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